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Abstract 

Background: Mouse is widely used in animal testing of cardiovascular disease. However, a large number of 
cardiovascular drugs that have been experimentally proved to work well on mouse were withdrawn because they 
caused adverse side effects in human. 

Methods: In this study, we investigate whether binding patterns of withdrawn cardiovascular drugs are conserved 
between mouse and human through computational dockings and molecular dynamic simulations. In addition, we 
also measured the level of conservation of gene expression patterns of the drug targets and their interacting 
partners by analyzing the microarray data. 

Results: The results show that target proteins of withdrawn cardiovascular drugs are functionally conserved 
between human and mouse. However, all the binding patterns of withdrawn drugs we retrieved show striking 
difference due to sequence divergence in drug-binding pocket, mainly through loss or gain of hydrogen bond 
donors and distinct drug-binding pockets. The binding affinities of withdrawn drugs to their receptors tend to be 
reduced from mouse to human. In contrast, the FDA-approved and best-selling drugs are little affected. 

Conclusions: Our analysis suggests that sequence divergence in drug-binding pocket may be a reasonable 
explanation for the discrepancy of drug effects between animal models and human. 

Keywords: Withdrawn cardiovascular drugs. Animal modeling. Sequence divergence. Side effects. Drug-binding 
pocket 



Introduction 

Mouse is the most commonly used vertebrate species in 
animal testing which is a vital part of drug development. 
Currently, all new pharmaceuticals undergo rigorous ani- 
mal testing before being licensed for human use. It is 
commonly accepted that animal testing on mouse can 
provide us critical information for assessing hazard and 
risk potential [1,2] . However, many withdrawn pharma- 
ceuticals worked very well in animal models while led to 
adverse side effects or failed to reach the expected effects 
in human. Over sixty drugs ever approved by the FDA 
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were withdrawn in the past twenty years [3]. A recent 
and well-known example is the failure of an Alzheimer's 
drug Dimebon (latrepirdine). The drug improved the 
performance of memory-impaired mice and rats [4,5], 
but failed to show a significant effect in the phase III 
clinical trial. Inevitably, a question arises as to why the 
discrepancy of drug effects between animal models and 
human. Answering this question will help us avoid drug 
withdrawals as much as possible and rescue' some of the 
withdrawn drugs by promoting genotype-based 
prescribing. 

We presumed that two aspects would be the prime 
determinants of the discrepancy of drug response be- 
tween animal models and humans. First of all, target 
proteins of withdrawn drugs may be functionally diver- 
gent between human and mouse. Phenotypic differences 
between species might be caused by changes in species- 
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specific interactions or gene expression [6-8], which also 
lead to the differential responses to the exogenous drugs. 
Secondly, three-dimensional structures of drug targets 
may have changed due to sequence divergence in the 
drug-binding pockets. As a result, although the drug tar- 
gets are functionally conserved, the binding patterns of 
drugs with these targets may change from animal models 
to human. 

In this study, both genetic aspects were examined to 
elucidate the mechanism of discrepancy of withdrawn 
drugs between human and mouse. We focused on with- 
drawn cardiovascular drugs, not only because cardiovas- 
cular drugs are in a position of importance in all 
pharmaceuticals, but also because the relevant informa- 
tion is comprehensive, including target, pharmaco- 
dynamics and so on. The results show that the drug 
targets are not involved in species-specific interactions. 
The binding patterns of withdrawn cardiovascular drugs 
are indeed affected by sequence divergence in the drug- 
binding pockets. The trend of binding affinity of these 
withdrawn drugs is to be reduced from mouse to human. 
Our study gleans valuable insights into the development 
of cardiovascular drugs and may reduce the amount of 
laboratory effort required in this process. 

Results 

Target proteins of withdrawn cardiovascular drugs are 
functionally conserved between human and mouse 

In order to compare the target protein of withdrawn car- 
diovascular drugs between human and mouse, we first 
examined the orthologous relationships of these proteins 
through Ensembl gene orthology prediction pipeline 
(http://www.ensembl.org/). It shows that all of the ortho- 
log relationships (25/25) are one to one (Additional file 
1: Table SI), which are widely assumed to have similar 
functions and cause similar phenotypes [9]. 

We further detected whether targets of withdrawn car- 
diovascular drugs were involved in species-specific inter- 
actions by reconstructing interaction networks for 
human and mouse drug targets separately through inte- 
grating experimental datasets (Materials and methods). 
To our surprise, there is only a narrow overlap between 
the reconstructed human and mouse drug target inter- 
action networks (Figure 1). Besides the drug targets, only 
four proteins (Nomenclature gene names are CCPllO, 
HDAC3, KCNMAl and PPARGCIA) are found in both 
networks. We supposed that it was mainly caused by in- 
adequacy of information in databases and confirmed 
each interaction by computational methods, such as lit- 
erature search, gene co-expression and bayesian net- 
works (Additional file 1: Table S2) [10]. The results 
indicate that most of the interactions related to with- 
drawn cardiovascular drug targets are highly preserved 
between human and mouse. 



In addition, we reasoned that changes in gene expres- 
sion underlie many of the phenotypic differences be- 
tween species and lead to differential responses to 
withdrawn drugs. As a result, we also measured the level 
of conservation of gene expression patterns of the drug 
targets and their interacting partners (Materials and 
methods). For the expression profiles of these ortholo- 
gous genes between mouse and human, are the average 
Pearsons correlation coefficient r is 0.27 and is signifi- 
cantly higher than that for the random gene pairs (Stu- 
dents two sample ^-test; p <10~^). It shows that the gene 
expression evolution of these drug targets between 
human and mouse was strongly shaped by purifying se- 
lection, suggesting that gene expression variation may 
not have an influence on the drug effect. 

Sequence divergence near the functional sites 

We reasoned that non-conserved substitutions of critical 
residues in drug targets could lead to different response 
to drugs. As a result, whether these residues had chan- 
ged from mouse to human was analyzed. The results 
shows that the critical sites of drug targets, such as active 
sites and phosphorylated residues, are conserved while 
residues nearby have changed a lot from mouse to 
human, including losses of hydrogen bond donors and 
substitutions of amino acid residues of opposite charge 
(Additional file 1: SI Table SI). For example, human 
plasminogen activator inhibitor 1, the target of troglita- 
zone (Drug Bank ID: DB00197), shares only 78% se- 
quence identity with its mouse ortholog. The critical 
sites are Arg369 and Met370 [11], which are conserved, 
while the residues within 10 A distance in three-dimen- 
sional structure have changed a lot, such as LysSllGly, 
His339Ser, and Glu373Thr. Previous work revealed that 
drug-binding sites on proteins usually exist out of func- 
tional necessity [12], so we can infer that the binding 
pockets of withdrawn drugs could be affected by se- 
quence divergence. 

Structural modeling of mouse drug targets and 
refinement by molecular dynamics simulation 

Because the structures of mouse drug targets are un- 
solved, we generated mutant models of human target 
proteins as rational mouse three-dimensional structures 
and employed molecular dynamics simulations for struc- 
tural refinement. According to previous work that MD 
simulations on the nanosecond time scale were sufficient 
for refinement of protein models [13], simulations for 2 
nanoseconds were first performed. The stability of each 
simulation was monitored through examination of struc- 
tural properties which occurred during the course of the 
2 nanoseconds production trajectories. The plots of 
root-mean-square-deviation (RMSD) from the original 
starting coordinates indicate that seven models are well- 
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Figure 1 Protein interaction networlcs of targets of withdrawn cardiovascular drugs. The interaction networl<s of human (panel A) and 
mouse (panel B) drug targets are displayed respectively. The pink nodes represent the drug targets, and the blue ones are proteins that have 
been experimentally proved to interact with drug targets. The edges between nodes represent interactions. The ellipse represents overlap 
between two networks. There are only four nodes to be found. 



behaved which means the simulations are reasonably 
converged and protein structures are well refined 
(Figure 2 A, C-E, G-I). Interestingly, a relatively large 
shift in RMSD (Figure 2F) at around 0.5 ns suggests that 
this structure is not well refined. It may be caused by the 
large loop (residues 1769-1787). The simulations of 
other two structures (Figure 2B, J) are also not well con- 
verged. As a result, we further performed simulations of 
these two structures for 5 nanoseconds until the RMSD 
values show only minor variation. 

The root mean-square fluctuations (RMSF) of the 
mutated residues in ten drug targets were also computed 
to examine the impact of mutating these residues on 
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Figure 2 RMSD values of mutant human protein structures 
during the dynamics simulations. The figures A-J respectively 
correspond to mutants of crystal structures (PDB code 1C5G, IDQA, 
IHWL, 1Y9C, 2BVR, 2KBI, 2QT9, 3D24, 3DZY and 3G43). 



conformational changes. The average RMSF of the 
mutated residues is significant higher than that of the 
conserved residues in 2 ns (Students two-sample ^-test, 
p = 3.52 X 10~^), indicating that mutated residues have a 
great influence on conformational changes of drug 
targets. 



Differences in the binding patterns of withdrawn 
cardiovascular drugs to their targets between human and 
mouse 

We next analyzed whether the drugs might bind human 
and mouse drug targets with different patterns and affin- 
ities. These drugs belong to different categories (Add- 
itional file 1: Table SI) and so are demonstrated 
separately. 

Troglitazone works as vasodilator agents and platelet 
aggregation inhibitors [14]. It selectively target six pro- 
teins, of which four have crystal structures (Nomencla- 
ture gene names: ESRRA, ESRRG, SERPINEl and 
PPARG). The interaction between troglitazone and 
Estrogen-related receptor gamma (ERR gamma-2) shows 
no interspecies difference because the protein sequences 
in human and mouse are identical. However, the other 
three targets that are not paralogs and perform different 
functions show variance in the drug-binding pockets. A 
few important amino acids within the binding pocket in 
human Peroxisome proliferator-activated receptor 
gamma (PPAR-gamma) have changed from mouse to 
human (Asn302Ser, Asn355Ser and His454Gln). It indi- 
cates that the troglitazone-binding pocket in mouse 
PPAR-gamma is distinct from that in human ortholog 
(Figure 3, A-B). The binding free energies demonstrate 
that the binding affinity of PPAR-gamma for troglitazone 
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Figure 3 Binding patterns of troglitazone with its targets. Panel (A) illustrates interactions between troglitazone and human PPAR-gamma, 
while (B) represents interactions between drug and corresponding mouse target. Panel (C) illustrates interactions between troglitazone and 
human PAI, while (D) represents interactions between drug and corresponding mouse target. The residue substitutions in the binding pocket are 
labeled. 



from mouse to human is reduced (Additional file 1: 
Table S3 and Materials and methods; Students two-sam- 
ple ^-test, p = S3x 10~^). As to plasminogen activator in- 
hibitor 1 (PAI), it shows that the residues in the binding 
pocket have changed greatly including, inter alia, 
Arg27Pro, HisSOTyr, ThrSlVal, Ser63Ala and so on 
(Figure 3, C-D). The binding energies also show that the 
specificity of interaction between troglitazone and recep- 
tor PAI is significantly different (Students two-sample t- 
test, = 1.14 X 10"^^). Unlike interacting with other tar- 
gets, troglitazone binds to Steroid hormone receptor 
ERRl as an inverse agonist. [15] There are only three 
amino acid substitutions of similar properties near the 
pocket (Additional file 1: Figure SI, A-B), and the recep- 
tor binding affinity is little affected (Students two-sample 
^-test, p = 0.078). 

Cerivastatin acts as antagonist on HMG-CoA reduc- 
tase, which functions as a homodimer and is the rate- 
limiting enzyme of sterol biosynthesis [16,17]. The com- 
putational docking results indicate that cerivastatin binds 
to the homodimer interface to keep it from dimerizing. 



The analysis shows that several critical residue substitu- 
tions (A: TyrSSlPhe, A: Glu548Gly and A: Ser574Gly), 
which lead to the loss of hydrogen bond donors, account 
for the difference in the binding pattern of cerivastatin 
between human and mouse (Figure 4, A-B). The analysis 
of binding free energy also indicates that the binding af- 
finity is reduced from mouse to human (Students two- 
sample ^-test, p = 0.033). 

Encainide can mediate the voltage-dependent sodium 
ion permeability of excitable membrane by blocking so- 
dium channel protein type 5 subunit alpha (HHl) [18]. 
Two residue substitutions (Leul812Ser and Alal813Val) 
occurred in the Encainide binding pocket (Additional file 
1: Figure S2, A-B). The amino acid residue Ser 1812 acts 
as a hydrogen donor according to the identification re- 
sult of drug-binding sites. However, the mouse and 
human HHl show similar encainide binding affinity. 
(Student s two-sample ^-test, p = 0.64). 

Mibefradil belongs to calcium channel blockers. Its tar- 
gets include T-type and L-type calcium channels, both of 
which mediate the entry of calcium ions into excitable 
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Figure 4 Binding patterns of cerivastatin with its targets. Panel (A) illustrates interactions between cerivastatin and human HMG-CoA 
reductase, while (B) represents interactions between drug and corresponding mouse target. The residue substitutions in the binding pocket and 
chains are labeled. 



cells and are also involved in a variety of calcium- 
dependent processes, such as muscle contraction, hor- 
mone or neurotransmitter release and so on [19,20]. Of 
the 11 mibefradil targets, only four proteins have peptide 
fragments resolved and one meets our standards (PDB 
ID code 3 G43). The computational docking results show 
that the binding pocket is composed of two helixes sep- 
arately from two chains, in which two different amino 
acids (Asnl607Gly, Alal608Ser) between species change 
the helix-helix interaction and so change the mibefradil 
binding pattern (Figure 5, A-B). The binding energy cal- 
culation also reveals that the binding affinity is reduced 
from mouse to human (Students two-sample ^-test, 
p = 1.53 X 10'^). 

Ximelagatran was investigated as an anticoagulant, 
which inhibited prothrombin from converting fibrino- 
gen to fibrin and activates factors [21]. The sequence 
identity of prothrombin between mouse and human 
was just 0.81, hinting that the drug-receptor interaction 




Figure 5 Binding patterns of mibefradil with its targets. Panel 
(A) illustrates interactions between mibefradil and human calcium 
channel, while (B) represents interactions between drug and 
corresponding mouse target. The residue substitutions in the 
binding pocket and chains are labeled. 



might also have interspecies specificity. Our results 
demonstrated that the residue substitutions (Glul49lys, 
Alal70Asp, Phel84Tyr, Lys222Asp and so on) in the 
binding pocket (Figure 6, A-B) changed the drug bind- 
ing pattern. The binding affinity of ximelagatran was 
also affected between species (Students two-sample 
^-test, p = 3,8x 10"^). 

The binding pattern and binding affinity is similar for 
FDA-approved drugs 

We next inspected that whether sequence divergence 
affects the binding pattern and binding affinity of FDA- 
approved drugs. In comparison, the FDA-approved and 
best-selling drugs, atorvastatin and clopidogrel (brand 
names are lipitor and plavix separately) were estimated. 
The results show that the binding patterns and binding 
affinities of both drugs are little affected by sequence di- 
vergence between human and mouse (SI Text Results). 

Discussion 

In this study, we demonstrate a reasonable possibility of 
the discrepancy of withdrawn cardiovascular drugs be- 
tween mouse model and human. The drug targets are 
proved to be not involved in species-specific interactions. 
Then, reasonable models for mouse protein structures 
are generated through molecular dynamics simulations. 
After that, we determine the differences of binding pat- 
terns and affinities of the withdrawn cardiovascular 
drugs between human and mouse and two FDA- 
approved and best-selling drugs are selected as controls. 
Our results show that the binding patterns of withdrawn 
cardiovascular drugs in our study are indeed affected by 
sequence divergence, especially the non-conserved resi- 
due substitutions, such as Lys into Glu, Ser into Ala and 
so on. The trend of binding affinity of these withdrawn 
drugs is to be reduced from mouse to human, which 
may explain the low specificity of human targets. Finally, 
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we explore whether off-target effects could be caused. It 
shows that more off-target effects can occur in human 
than in mouse due to the local structural similarities of 
drug-binding sites. 

Usually, adverse side effects occur when the drug inter- 
acts with unintended targets [22]. Since the binding sites 
and affinity of withdrawn cardiovascular drugs were esti- 
mated to have changed from mouse to human, these drugs 
might interact with more off-targets in human than in 
mouse, due to similar local structures of drug-binding 
sites. So, we determined the cardiovascular drug-binding 
sites of human and mouse targets separately and then 
searched the similar local structures in structure database. 
It demonstrates that human targets of withdrawn drugs 
get 7 off-targets in all while mouse targets get no off-tar- 
gets due to the differences of drug-binding sites (Materials 
and methods; Additional file 1: Table S4). Adverse side- 
effects may arise due to significant sequence divergence 
between mouse and human (Figure 7). Before drug testing, 
the animal model should be selected in a careful and me- 
ticulous manner, not only in CAD systems but also in 
other complex disease systems. 

A key issue in translational medicine is the need to de- 
vote more attention to the characteristics of the drug- 
target interaction [23], yet most of target proteins in ani- 
mal models have not or fragmentally solved. Besides, 
according to previous works, one missense mutation is 
sufficient to lead to altered binding pattern of drug tar- 
gets [24,25]. Studies on structural divergence of drug tar- 
gets between disease models and human provide 
valuable information for the development of human dis- 
ease models. In this study, we employed MD simulations 
to construct models of good quality as target structures 
and then applied them to comparison of drug-target 
interactions. It will be very useful to identify a suitable 
model for animal testing, in case of false positive results 
supported by inappropriate experiments. 



Although we highlight the contributions of sequence 
divergence, we also illustrate specific interactions that 
were experimentally confirmed in just one species. Re- 
cent studies have shown that species -specific interactions 
and regulations exist in human and mouse [26,27]. These 
species-specific properties may lead to different 
responses to the same pharmaceutical. Besides, false 
positive results in mouse may also lead to the success of 
animal testing while failure of clinical treatments. For ex- 
ample, prenylamine is a calcium channel blocker of the 
amphetamine chemical class which is used as a vasodila- 
tor in the treatment of angina pectoris. It has been 
shown to partially metabolize to amphetamine and can 
cause false positives for it in drug tests [28-30]. 

In summary, our results provide a reasonable explan- 
ation for the discrepancy of cardiovascular drugs. We 
will further testify for other pharmaceuticals. 



Mouse 



Human 



Target Drug 




Off-target 



Off-target 



Figure 7 Illustration of side-effects caused by sequence 
divergence. Adverse side-effects may arise due to significant 
sequence divergence in tine drug-binding pocl<ets between mouse 
and liuman. 
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Materials and methods 

Collection of withdrawn and FDA-approved cardiovascular 
drugs 

The information of withdrawn drug for CADs were 
retrieved from DrugBank database [31] and previous works 
[3]. Seven withdrawn cardiovascular drugs (Drug Bank ID: 
DB00197, DB00439, DB01228, DB01388, DB04825, 
DB04831 and DB04898) were included. In this dataset, pre- 
nylamine (Drug Bank ID: DB04825) was withdrawn for it 
partially metabolized to amphetamine and could cause false 
positives for it in drug tests [29]. It was just a prodrug, so it 
was excluded out of our study. Ticrynafen (Drug Bank ID: 
DB04831) was also excluded because it had no target infor- 
mation. We also chose two FDA-approved and the most 
sales-potential drugs, Lipitor (Drug Bank ID: DB01076) and 
Plavix (Drug Bank ID: DB00758) as positive control because 
the higher sales-potential might indicate the more satisfac- 
tory efficacy and less side-effect. 

Reconstruction of interaction networks of withdrawn drug 
targets and literature search 

The interaction data was retrieved from several common 
databases, including DIP [32], BIND [33], HPRD [34], 
BioGRID [35], MINT [36] and IntAct [37]. A Cytoscape 
[38] plugin, BisoGenet, was applied to reconstruct inter- 
action networks of human and mouse drug targets separ- 
ately. Only the interactions that were verified by 
experiments were displayed. Another cytoscape plugin, 
Agilent Literature Search 2.71, was used to do literature 
search. Gene aliases and context were set to restrain the 
search results. 

Analysis of gene expression data 

The microarray datasets (GSE8000 for mouse [39] and 
GSE 11560 for human [40]) were downloaded from the 
Gene Expression Omnibus (GEO) (http://www.ncbi.nlm. 
nih.gov/geo). The expression profile of the drug targets 
were analyzed according to our previous work [41]. 

Preparation of structures of drug targets 

All the protein structure files of withdrawn drug targets 
were downloaded from PDB database [42]. Structures, 
containing the main functional domains and of which se- 
quence coverage is above 80%, were adopted to do the 
subsequent analysis. If the sequence coverage was below 
the cutoff value and yet the drug-binding regions were 
intact and clearly reported in a structure, then the struc- 
ture was also included in the dataset (Additional file 1: 
Table SI). 

Generation of mutant models for mouse targets and 
refinement 

Mutant models of human protein structures were made 
as mouse protein structures using MODELER program. 



[43] High levels of molecular dynamics with simulated 
annealing were performed to optimize the local struc- 
tures of mutant models. After that, molecular dynamics 
simulations were performed using NAMD and 
CHARMM31 force filed with CMAP correction [44,45]. 
The protein models were solvated in a TIP3P water box 
[46] and ionized by NaCl (0.152 M) to mimic physio- 
logical conditions. The ionized systems were minimized 
for 50,000 integration steps and equilibrated for 2 ns 
with 2 fs time stepping. Following this, a 2 ns uncon- 
strained equilibration was performed for subsequent tra- 
jectory analysis, with frames stored each picosecond. 
Constant temperature (T = 310 K) was enforced using 
Langevin dynamics with a damping coefficient of 5 ps. 
Constant pressure (p=l atm) was enforced through the 
Nose-Hoover Langevin piston method with a decay 
period of 100 fs and a damping time constant of 5 per 
picosecond. Van der Waals interaction cutoff distances 
were set at 12 A, (smooth switching function beginning 
at 10 A) and long-range electrostatic forces were com- 
puted using the particle-mesh Ewald (PME) with a grid 
size of less than 1.0 A. 

Comparison of binding patterns of cardiovascular drugs 

We compared modes of actions of withdrawn cardiovas- 
cular drugs following the steps below. Step one; we de- 
tect the drug binding pockets through a CHARMm- 
based molecular dynamics scheme [47]. We performed 
simulated annealing refinement and saved the top ten 
poses to prepare for the calculation of binding energy. 
Step two: a flexible docking method was applied to find 
the optimal binding sites by simulating protein flexibility 
and docking drugs with an induced fit receptor 
optimization [48]. Step three; we determined the critical 
drug-binding sites by using Receptor-Ligand Interactions 
tool. We set 2.3 A as cutoff value to exclude the minor 
important bonds. All simulations reported above were 
carried out using Discovery Studio© v2.5.0.9164, built on 
the SciTegic Enterprise Server platform (Accelrys Inc, 
2009). 

Detection of local structural similarities of drug-binding 
sites in structure databases 

An accurate algorithm for detecting local protein struc- 
tural similarity, CMASA [49] was applied to detect local 
structural similarities of drug-binding sites in non-redun- 
dant structure database from SCOP database [50]. The 
binding sites detected above were used as queries. Sub- 
stitution files were applied in the process and the cutoff 
was set to the default value. 

Binding free energy calculations and statistical analysis 

The drug-target binding free energies were calculated in 
the protocols of Discovery Studio, according to the 
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previous work [51]. The values of the binding free energy 
(AG binding) for each drug-target complex were calcu- 
lated based on the following equation: 

AG binding = AG complex - AG drug - AG target 

(1) 

We used the top ten poses mentioned above to calcu- 
late binding energies. In situ ligand minimization and lig- 
and conformational entropy were set to default value. 
The distance cutoff value of ligand conformational en- 
tropy was set to 14.0 A. If the ten calculated binding en- 
ergies (Additional file 1: Table S3) followed a normal 
distribution, Students two-sample t test was applied to 
compare the binding affinity of drugs to their targets be- 
tween human and mouse. 

Additional file 



Additional file 1: Supporting Information. 
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